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We report the direct detection of radiation pressure on the asteroid 2009 BD, one of the smallest 
multi-opposition near-Earth objects currently known, with H ~ 28.4. Under the purely gravitational 
model of NEODyS the object is currently considered a possible future impactor, with impact solutions 
starting in 2071. The detection of a radiation-related acceleration allows us to estimate an Area to Mass 
Ratio (AMR) for the object, that can be converted (under some assumptions) into a range of possible 
values for its average density. Our result AMR = (2.97 ± 0.33) x 10 -4 m 2 kg -1 is compatible with the 
object being of natural origin, and it is narrow enough to exclude a man-made nature. The possible origin 
of this object, its future observability, and the importance of radiation pressure in the impact monitoring 
process, are also discussed. 



1 Introduction 

The search for potentially hazardous asteroids 
has yielded the discovery of many objects in orbits 
very similar to the orbit of the Earth. Frequently 
such discoveries trigger the question as to whether 
the object is of natural or man-made origin, which 
is usually addressed by integrating the orbit back- 
ward in time to see if a close Earth approach oc- 
curred at a time consistent with the launch of a 
rocket into a heliocentric orbit. 

A variety of forces act on solar system objects, 
including gravity, outgassing, solar radiation pres- 
sure, the Yarkovsky effect, the Poynting- Robertson 
effect, even possibly electrostatic forces. For the 
vast majority of objects, their motions can be satis- 
factorily explained by considering only the gravita- 
tional force. Even though several of these forces can 
be acting on an object simultaneously, the contri- 
bution from the weaker forces may be undetectable, 
and therefore no attempt is made to include them 
when solving for the dynamics of the object. When 
the observed motion of an object departs from 
what one would expect from a purely gravitational 
model, additional forces need to be taken into con- 
sideration. Which of these forces gets taken into 
account can depend on the observed nature of the 
object. For example, it is natural to include the ef- 
fects of outgassing on an object showing cometary 
activity. On the other hand, a small apparently 
asteroidal object is a candidate for the inclusion 
of solar radiation pressure. In such cases, its de- 



tection can yield the bulk density of the object, 
thereby making it easy to distinguish between the 
two possible origin cases. One such object, known 
by the internal designation assigned by the initial 
observer, J002E3, has now been confidently linked 
to the Apollo 12 fourth staged. This paper dis- 
cusses another such object, 2009 BD. 



2 Observational History 

2009 BD was discovered by the Catalina Sky Sur- 
vey, Mt. Lemmon station, on 2009 January 16 UT 
[1]. Over the next ten days, the object was ob- 
served astrometrically 119 times from 15 different 
stations. A few days later, the apparent magnitude 
became too faint and the solar elongation too small 
for continued observation during that apparition. 

The derived visual absolute magnitude was 
reported as H = 28.8, which corresponds to a 
diameter of about 7 m if one assumes a visual 
geometric albedo of 0.12, a rough average for 
asteroids in the inner portion of the solar system. 
The resulting orbit showed a semimajor axis of 
1.01 ua, an eccentricity of 0.03, and an inclination 
with respect to the mean ecliptic of J2000 of 1.5°. 
Its low Av with respect to the Earth of 3.73 km/s 
places it at the top of the list of near-Earth 
asteroids accessible with a spacecraft fry-by^ 2 ). 



^'www. jpl .nasa. gov/releases/2002/release_2002_178 . cfm 
("'echo . jpl .nasa.gov/~lance/delta_v/delta_v. flyby.html 
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We recovered 2009 BD at Mauna Kea Observa- 
tory with two 300 s exposures taken on 2009 Octo- 
ber 18 UT using the University of Hawaii 2.24 m 
telescope, with a confirming observation on 2009 
November 1 UT. On the latter date, it was nec- 
essary to stack eighteen 60 s exposures to reach 
R = 22.7 due to the 98% illuminated Moon being 
just 12° away from the object; background levels 
reached 44000 DN in these unfiltered, two-by-two 
binned images. Only a single position from the 
stacked image was reported for this night. 

These original astrometric solutions were per- 
formed using the USNO-B1.0 reference catalog 
[2], with proper motion for each source taken into 
account. The inclusion of these new observations 
by the impact monitoring sites did not eliminate 
the possibility of a future Earth impact. As of 
this writing, NEODyS shows five possible impacts 
in the next century^ 3 ), while Sentry shows eight 
possible impacts^ 4 ), the earliest being in 2071 in 
both cases. 

We observed 2009 BD at a third opposition with 
two 150 s exposures taken on 2010 July 14 UT, 
again with the UH 2.24 m telescope on Mauna 
Kea. On 2010 September 1 UT we obtained an- 
other observation with the 3.6 m Canada-France- 
Hawaii Telescope on Mauna Kea; the object had 
magnitude V ~ 23.5, and we obtained a SNR of 
6.6 with a single 300 s exposure. 

However, a gravitational orbit solution that in- 
corporated these new observations revealed residu- 
als far in excess of our expected error. Such situa- 
tions are not uncommon in our experience, and we 
usually attribute the problem to astrometric refer- 
ence catalog biases [3]. However, after repeating 
the solutions using the 2MASS catalog [4], which 
is believed to have substantially smaller biases [5], 
the declination bias actually became worse. At 
this point we began investigating the importance 
of non-gravitational forces acting on the object. 

3 Methods 

Our initial orbital computation used 119 astro- 
metric observations collected by various stations in 
the few days after the discovery, from 2009 January 
16 to 26, plus 6 high-precision positions obtained 
by us from Mauna Kea in Hawaii, spanning from 

'newt on. dm.unipi . it/neodys/ index .php?pc=l . 1 . 2&n=2009S 
'neo . jpl . nasa . gov/risk/2009bd . html 



2009 October 18 to 2010 September 1. 

The first part of the analysis was performed us- 
ing the radiation pressure option in the software 
Find_Orb( 5 ); it allows the computation of an or- 
bit while adding a seventh free parameter to the 
usual six orbital elements, namely the Area to Mass 
Ratio (AMR) of the object. The force acting on 
the object because of the incoming solar radia- 
tion is directly proportional to its cross-sectional 
area exposed to the Sun, and its effect on the or- 
bit (that is, the resulting acceleration) is there- 
fore proportional to the ratio between this cross- 
section and the mass of the object, that is the 
AMR. We did not consider the effects of out- 
gassing because there is no obvious cometary ap- 
pearance in any of our images, the object is too 
large for the Poynting-Robertson effect to be sig- 
nificant, and electrostatic forces would only be im- 
portant if there were some mechanism to make the 
object significantly charged. Because of the small 
size and likely rapid rotation [6] of 2009 BD, the 
surface of the asteroid may be sufficiently isother- 
mal to render the Yarkovsky effect undetectable. 

The AMR can then be combined with an esti- 
mate of the diameter of the object, obtained from 
its absolute magnitude H and from assumptions 
on its geometric albedo, to get a range a possible 
densities. A density value around or above that ex- 
pected for a porous rock (p ~ 1000 kg/m 3 ) would 
suggest a natural origin for the object, while a sig- 
nificantly lower value (p ~ 20 kg/m 3 ) would point 
toward a man-made origin (a hollow metallic shell, 
such as an upper stage of a rocket). 

Both parts of the process (AMR and diame- 
ter determination) are based on the assumption 
that the body is approximately spherical. The im- 
plication of a non-spherical shape, together with 
lightcurve constraints on its importance, are briefly 
discussed at the end of this work. 

3.1 Outlier rejection techniques 

The computed value for the AMR can be sen- 
sitive to each observation included in the solution, 
and even a few biased astrometric data points may 
be enough to skew the results, leading to an in- 
correct evaluation of the importance of radiation 
pressure on the best-fitting orbit. 

At the same time, the diameter computation re- 
lies heavily on the absolute magnitude H; however, 
no calibrated photometry of the object is available, 

'www. pro j ectplut o . com/f ind_ orb .htm 
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and it is therefore necessary to extract our best es- 
timate of H from the sparse photometry reported 
by various stations together with their astrometric 
positions. 

Since both of these results are sensitive to the 
input data, it is important to ensure that the 
dataset is internally consistent, and all outliers 
have been rationally rejected, both from the as- 
trometric positions and from the magnitude values. 

A commonly used method for outlier rejection, 
known as Chauvenet's criterion [7], may not appli- 
cable in this case, since it is possible that more than 
one datapoint will need to be rejected, while the 
criterion is statistically founded only when testing 
the rejection of a single outlier. Fortunately, an- 
other less known rejection method, named Peirce's 
criterion and first introduced in 1852 by Benjamin 
Peirce [8], is statistically correct even in case of 
multiple rejections (at the cost of a slightly higher 
mathematical complexity) . 

For a simplified (but mathematically equivalent) 
description of the method we refer the reader to 
[9] , which gives a more astronomy-oriented presen- 
tation of the criterion, together with formulas to 
compute its rejection thresholds in practical cases 
with a different number of free parameters. Un- 
fortunately the formulas, summarized in App. A, 
are not analytical, and they require an iterative 
process to obtain a numerical value. However, they 
can easily be implemented in a computer code, 
that can give sufficiently accurate values with only 
a few iterations. 

The criterion, in its original form, only assumes 
that the data points are normally distributed 
around their means, with the same standard devi- 
atW 6 ). We can therefore directly apply Peirce's 
criterion to our magnitude rejection problem, since 
it is safe to assume that the residuals of each 
magnitude are randomly distributed around the 
mean, and they reasonably follow a unidimensional 
Gaussian distribution. 



The last requirement is obviously not satisfied in this case, 
since every station produces astrometry and photometry 
with different accuracies. However, we can estimate their 
typical error bars from other astrometric and photometric 
data they reported, and then divide the residuals by this 
expected accuracy; every resulting residual will then be adi- 
mensional, and distributed as the standard N(0, 1) Gaussian 
that we need for the rejection criteria. This is equivalent to 
saying that we will perform our rejection tests in terms of 
number of sigmas, not in terms of absolute errors. 



Unfortunately, the same is not true for the re- 
jection of astrometric observations; in this case 
the meaningful metric for accuracy is the angu- 
lar distance between the observed position and the 
ephemeris from the fitted orbit, that is the RMS 
of the right ascension and declination residuals^ 7 ). 
This quantity, however, is positive definite, and it 
is obviously not normally distributed around zero; 
under the assumption that each component of the 
residual (right ascension and declination) is nor- 
mally distributed, the RMS is known to follow the 
Rayleigh distribution (a sort of square root of the 
chi-square distribution with two degrees of free- 
dom). This distribution can be used to compute 
a new version of Peirce's criterion that is valid in 
our case. The mathematical derivation of this re- 
sult is explained in App. A. 

3.2 Photometric debiasing and outlier 
rejection 

Most of the physical results obtained in this work 
are sensitive to the assumed value for the abso- 
lute magnitude H of 2009 BD. Unfortunately, no 
accurately calibrated photometry of this object is 
known to exist, and we are forced to use the poorly 
calibrated and temporally sparse photometry as- 
sociated with each astrometric position submitted 
to the Minor Planet Center (MPC) by various sta- 
tions. 

Such photometric values are known to be affected 
by various sources of error, both statistical and sys- 
tematic. It is therefore necessary to perform an 
accurate analysis of each one of them, to properly 
estimate its bias with respect to a hypothetical cal- 
ibrated photometry, its weight in the overall av- 
erage and the possibility of its complete rejection 
from the analysis. The complete procedure used 
to compute our best value of H (in the usual H- 
G formalism from [10]) is quite complex, and it is 
described in detail in App. B. 

The application of Peirce's criterion in this case 
leads to the rejection of photometric data from 3 
nights of the discovery opposition; two of those 
rejected points, both from station J95, were taken 
at reasonably high phase angles (a = 52° and 

^ 7 'We preferred not to use the two coordinates as independent 
quantities because we wanted to avoid the need of rejecting 
only one coordinate, and keeping the other. Each astromet- 
ric observation is treated as a single entity, and it is either 
rejected or accepted entirely. 
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to extract high-precision astrometry from trailed 
detections. A single 2MASS-based position was 
obtained from a stack of all Tzec Maun exposures, 
with a total SNR around 25. We replaced all the 
other MPC-reported astrometry from code H10 
with this single position (also listed in Tab. 1) and 
used it in our computations. 
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a = 77°); this geometry may explain the anoma- 
lous offset, either because of a bad assumption 
of G = 0.15, or because of the more extreme 
light curve amplitude typical of high phase angle 
observations [11]. 

The remaining photometric values are then used 
to compute our best value for H, and its error bar. 
We obtain a value of H = 28.43 ± 0.12, where the 
quoted uncertainty is purely statistical, and does 
not take into consideration the lightcurve ampli- 
tude. 

3.3 Astrometric outlier rejection 

As explained above, the results of this work 
are also extremely sensitive to the astrometric 
dataset used as a starting point in the orbital 
computation, and for this reason we must be sure 
that all the astrometry is cleaned from biases and 
systematic effects that can skew the final results. 
To attain this result, we tried to both maximize 
the accuracy of the astrometric positions that 
we produced, and at the same time we used the 
methods discussed above to reject suspect data 
points from the astrometry of other stations on 
which we don't have direct control on the data 
reduction process. 

First of all, our astrometric positions from 
Mauna Kea were all re-reduced against the 2MASS 
catalog, which is currently considered the least 
affected by systematic biases in the magnitude 
range we need. Since 2MASS does not have 
magnitude information in the visible range, the 
magnitudes from the USNO-B1.0 reductions were 
kept and used in the photometric analysis above. 
Our astrometric and photometric observations are 
listed in Tab. 1. 

We also noticed that a set of eight positions from 
Tzec Maun Observatory (observatory code H10) 
showed unusually large residuals; we found that 
the FITS files of those observations were placed 
online by the observers, and a total of ten indi- 
vidual exposures were available for download®; 
since the object appeared significantly trailed and 
faint in all those frames, we re- measured them ( 9 ) Since a variety of different catalogs were used by the various 
with our astrometry tools, specifically designed observatories, and the positions span a range of more than 

70° in the sky, we will focus only on rejecting outliers, not 

on locally debiasing each specific astrometric position. 



Table 1 : Astrometric position and photometry for 
our 6 observations, plus the stacked re- 
reduction from the H10 data. The as- 
trometry is referred to 2MASS, while 
the magnitudes are calibrated against 
USNO-B1.0. 

The second part of our strategy targeted the ob- 
servations obtained by other stations, all grouped 
during the first 10 days after the discovery in 2009. 
Being clustered in a short timespan, they carry less 
weight compared to ours in the final solution, but 
it is nevertheless necessary to ensure that at least 
their average is not biased in any particular direc- 
tion, which may mimic an offset introduced by ra- 
diation pressure effects®. 

We applied our "Rayleigh version" of Peirce's cri- 
terion (discussed in App. A) to the full astromet- 
ric dataset, including our observations from Mauna 
Kea. 

First of all, for each observation we computed an 
expected astrometric error bar using the historical 
performance of the site as published by the Minor 
Planet Center ( 10 ). The two-component residuals 
published by the MPC were combined together in 
quadrature, to obtain a single error estimate that 
we assume to have Rayleigh distribution. For our 
Mauna Kea observations, we used an a-priori error 
of 0.13", an average of the error bars derived during 
our astrometric solutions^ 11 ). For the remeasured 



f its . tzecmaun2 . org/temp/2009bd_f its . zip 



(8) From 

publicly posted by R. Wodasky on mpml 
tech. groups . yahoo . com/group/mpml/message/21487 



at ' 10 ' www. minorplanet center . org/ iau/ special/residuals . txt 

' 11 'The error bars derived by the astrometric software range 
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position from station H10 we assumed an error of 
0.25", from our astrometric software. 

This expected accuracy for each station is used 
twice in the overall orbital computation. First of 
all, it is converted into a weight in the software 
Find_Orb, to properly account for each station's 
performance in the overall orbital fit (in the absence 
of any better information on the astrometric error). 

We then obtain a first-pass orbit using all 
available observations, with the seven-parameter 
radiation-pressure model. This orbit gives us resid- 
uals for each astrometric position, and they are 
then normalized with the error value computed 
above for each station, as required by Peirce's 
statistic. 

At this step Peirce's criterion is applied (in 
its Rayleigh form) to reject all possible outliers 
from the astrometry set. A rejection threshold is 
computed for the case of N = 118 (12) datapoints 
and m = 7 fitted parameters, and for various 
numbers of possible rejections n, using the recur- 
sive relations in App. A. Each value of n is then 
tested against the data, finding the first value n 
where less than n observations have normalized 
residuals greater than the corresponding computed 
threshold. In our sample that happens at n = 11, 
meaning that the 10 observations with largest 
normalized residuals need to be rejected from the 
final orbital solution. 

This new dataset, now including only 108 astro- 
metric positions (our 6 from Mauna Kea, that were 
all kept by the criterion, plus the 102 from other 
stations that survived the rejection process), will 
be the basis from which our final AMR value will 
be computed. 

4 Results 

This dataset of 108 astrometric positions can 
now be fitted with the software Find_Orb, us- 
ing the same weighting scheme described above. 
The best-fitting purely gravitational solution has 
an RMS of 0.395" and a reduced chi-square of 

between 0.11" and 0.15", and are mostly dominated by an 
assumed 0.10" catalog bias. This assumption may be conser- 
vative for the 2MASS catalog, but we wanted to make sure 
that our results were not affected by an excessive weighting 
of our observations. 

The total number of observations went from TV = 125 to 
N = 118 because of the replacement of 8 positions from 
H10 with our single stack. 



Xu = 1.390, and it shows unreasonably large resid- 
uals for some of the Mauna Kea observations (see 
Tab. 2). 

We repeated the orbit solution allowing for so- 
lar radiation pressure but otherwise leaving the as- 
trometric data and weights unaltered. The result- 
ing RMS residual was only 0.291", 71% as large 
as for the gravitation-only orbit, with a significant 
improvement in the residuals. The reduced chi- 
square went down by a factor of more than two to 
Xu = 0.626, indicating both the significance of the 
additional parameter and that we were conserva- 
tive in our estimates of the astrometric uncertain- 
ties. The best-fit area to mass ratio associated with 
this solution is 2.97 x 10~ 4 m 2 kg -1 . 

The fact that the xl statistic is not signifi- 
cantly greater than unity indicates that our radia- 
tion pressure model is sufficient to satisfy the ob- 
servations. We therefore do not need to include 
the Yarkovsky effect, but that does not mean that 
the effect is not present, just that it remains unde- 
tectable in the current dataset. 



Date 
[UT] 


A«["] AS["] 
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Aa["] AS["} 
Non-grav. 


2009-10-18.429045 
2009-10-18.433510 


-0.01 -0.07 
-0.02 -0.12 


-0.03 -0.04 
-0.05 -0.09 


2009-11-01.468038 


-0.19 +0.09 


-0.09 +0.17 


2010-07-14.587053 
2010-07-14.589045 


+0.26 -0.09 
+0.28 -0.09 


+0.09 -0.01 
+0.11 -0.01 


2010-09-01.322573 


-0.53 -0.13 


-0.14 -0.06 



Table 2: Comparison of the residuals of the purely 
gravitational solution with those includ- 
ing the radiation pressure effect. Some 
of the gravitational residuals, such as 
a 0.53" residual for the 2010 Septem- 
ber 1 observation, are many sigmas away 
for our best-estimate for their error bars 
(typically around 0.13" or less for a 
2MASS reduction). 

We then generated 10 000 sets of synthetic 
astrometric observations by adding Gaussian noise 
to the 108 available observations. The noise was 
generated assuming again a standard deviation 
equal to the RMS of the station, computed from 
the MPC residuals as explained above. Orbit 
solutions allowing for solar radiation pressure were 
performed on each set, yielding an uncertainty in 
the area to mass ratio of 0.33 x 10~ 4 m 2 kg -1 , 
which represents a 9a detection^ 13 ) of that non- 

The significance quoted here is purely statistical, and it does 
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gravitational effect. 

We can now use this estimate of the object's 
AMR = (2.97 ± 0.33) x 1(T 4 m 2 kg" 1 , together 
with our best estimate of its absolute magnitude 
obtained above (H = 28.43 ± 0.12), to put some 
constraints on the object's density. 

If we first assume an albedo of 0.12, the typi- 
cal average for asteroids in the inner portion of the 
solar system, we get a density of p ~ 640 kg/m 3 , 
which is substantially higher than the densities of 
the two man-made objects for which similar anal- 
yses have been performed and consistent with the 
density of very porous rock. The Apollo 12 fourth 
stage (essentially an empty fuel tank) has a den- 
sity around 20 kg/m 3 , while the density of the ob- 
ject designated as 6Q0B44E was measured at about 
15 kg/m 3 ( 14 ). 

In comparison, the NEAR-Shoemaker spacecraft 
flyby of (153) Mathilde yielded a bulk density 
determination of 1300 kg/m 3 . Although such 
densities are less than what would be expected 
for solid rock, the fact that Mathilde did not 
shatter despite the impacts that left craters with 
diameters a significant fraction of the size of the 
asteroid strongly suggests an unconsolidated inte- 
rior structure, with porosity approaching 50% [12]. 
Furthermore, the object 6R10DB9, now recognized 
as natural and designated 2006 RH120, has a 
measured density around 400 kg/m 3 , comparable 
with our value for 2009 BD. 



4.1 Effect of the unknown albedo 

Of course, the density of 2009 BD depends 
strongly on the visual geometric albedo, for which 
we have no measurement. It may therefore be 
more meaningful to explore the density values 
associated with a broad range of albedoes. This 
analysis is summarized in Fig. 1, showing that the 
best-fit density stays well above the typical value 
for a man-made object down to albedoes smaller 
than 0.01. 

Although the asteroid was near 90° solar elon- 
gation at about the time that WISE began its all- 
sky survey in 2010 January, and then again in late 

not take into account possible systematic effects, such as the 
omission of some other non-gravitational effect. 
Computed from the AMR estimate available at 
home .gwi .net/~pluto/mpecs/6q0b44e .htm, assuming 
an albedo similar to that used for 2009 BD. 



Density [kg /in 3 ] 

2000 - 




Figure 1: Estimate of 2009 BD density for albedo 
values ranging from extremely dark 
(pv = 0.01) to almost complete reflec- 
tion (p v = 1.00). 

June of the same year, the object was not detected. 
Ephemeris uncertainties were less than 2" at the la 
level, so the fact that no WISE source was associ- 
ated with a known solar system object cannot be 
attributed to large ephemeris error. 



5 Conclusions 

The astrometric data reported in this work show 
the necessity of modeling the orbit of 2009 BD with 
a perturbation caused by a non-gravitational effect, 
capable of offsetting the sky-plane position of the 
object of at least 0.5" during an observed arc of less 
than 2 years. 

We propose that this offset can be explained 
with the addition of a single parameter, namely 
an AMR = (2.97 ± 0.33) x 10~ 4 m 2 kg" 1 , caused 
by the action of radiation pressure on the ob- 
ject. This result is in agreement with the value 
AMR ~ 1.9 x 10" 4 m 2 kg~* that would be ex- 
pected for a spherical rocky object of this absolute 
magnitude (H = 28.43±0.12), assuming water den- 
sity and an albedo of 0.12. 

This agreement is consistent with the hypoth- 
esis that the effect we are detecting is related to 
the radiation pressure acceleration. Other non- 
gravitational effects may contribute on a smaller 
level. In particular, the multi-parameter Yarkovsky 
effect is not needed to explain the current dataset. 
The Poynting-Robertson effect is also negligible for 
an object of this size, being a factor of ~ ^ ~ 10~ 4 
smaller than the radiation pressure effect. 
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5.1 Possible origin 

As already mentioned, the density range derived 
above points toward a natural origin for 2009 BD, 
excluding the possibility of a man-made object, 
such as rocket debris, that was left in a heliocentric 
orbit at the time of the launch. It is therefore 
necessary to explore other scenarios that may be 
capable of placing a natural rocky object of this 
size into an Earth-like orbit. 

The most direct explanation would be that 
2009 BD is a regular NEO coming from the Main 
Belt through the usual v§ resonance mechanism; 
however, it is unlikely that such a mechanism 
would produce an orbit resembling that of the 
Earth, except just by coincidence. Some works, 
such as [13], suggest that there may be an enhance- 
ment of objects delivered to a low-eccentricity 
Earth-like orbit with this mechanism, but the 
overall likelihood of an object being discovered in 
that state is low (and even lower if we consider that 
many other objects, such as 1991 VG, 2000 SG344 
or the already mentioned 2006 RH120, also share 
similar Earth-like orbital properties). 

More likely, this object could be a fragment of 
the Moon ejected during a collision and cratering 
process on our satellite. This theory, already pro- 
posed to explain the Earth-like orbit of other ob- 
jects, starting with 1991 VG [14], can directly ex- 
plain most of the orbital properties of 2009 BD. 

In particular, fragments ejected from the Moon 
with the lowest possible escape velocity (v ~ 
2.4 km/s) can enter heliocentric orbits passing 
through the Earth Lagrangian points; in such cases, 
the resulting orbit [15] closely resembles that of 
2009 BD, with a semimajor axis slightly larger 
or smaller than the Earth's, eccentricities < 0.1 
and inclinations < 1°. The current inclination of 
2009 BD may look larger than this value, but back- 
ward integrations show that its eccentricity was 
well below that threshold until a close encounter 
with the Earth in 1955. 

The main caveat of this hypothesis is that the 
size of 2009 BD is close to the upper limit of what 
can probably be lifted (and achieve escape veloc- 
ity) from the Moon surface during an impact event. 
However, studies [14] for the case of 1991 VG (an 
object of similar size, H ~ 28.4) confirm that 
the scenario is plausible and reasonably likely on 
timescales comparable with the dynamical survival 



time of an object on such an orbit. 

The Moon ejection scenario may be confirmed 
or disproved with albedo or spectral information, 
that will hopefully be collected during the next 
apparition of 2009 BD. 

Finally, it may be worth mentioning at least two 
more scenarios, probably much less likely than the 
one discussed above. 

A direct generalization of the lunar ejection case 
would be to assume that 2009 BD is a fragment 
of the Earth surface ejected during a past collision; 
this possibility is unlikely because the large ejection 
velocities required to achieve an heliocentric orbit 
from the Earth's surface would require a massive 
impact on a timescale less than 10 6 yr ago (the av- 
erage survival time of objects with this kind of or- 
bits), that would have likely left obvious geological 
signs on our planet. 

A second hypothesis, that may be worth explor- 
ing if the Moon ejection scenario becomes unlikely 
due to physical observations, is the possibility that 
the object entered this Earth-like orbit after a very 
close encounter with our planet, by slowing down 
significantly by friction during a fly-by phase inside 
the Earth atmosphere. Similar events have been 
observed from the ground, and they may lead to 
Earth-like orbits such as that of 2009 BD. 

5.2 Implications for Impact Monitoring 

2009 BD is currently classified as a possible 
impactor by both NEODyS and JPL's Sentry, 
with low-probability impact solutions starting from 
2071. The small size of the object, around 8 m 
assuming an albedo of 0.12, makes it a negligible 
threat to the Earth, being at most capable of caus- 
ing an event similar to the 2008 TC3 impact, with 
almost no ground damage and only a few meteoritic 
fragments reaching the ground. The situation can 
be slightly different if the object is of metallic com- 
position; in that case the bulk of the mass may 
be able to reach the ground, but it would still be 
too small to cause any significant and widespread 
threat or damage. Furthermore, our derived bulk 
density value clearly points toward a very porous 
rock, and it is incompatible with a metallic com- 
position, making the scenario of a ground impact 
unlikely. 

However, it is still interesting to explore the im- 
plication of the radiation pressure effect on the 
orbital evolution of the object. The best-fitting 
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gravitational and non-gravitational solutions are 
now only about 20 km apart, but they will tend 
to diverge in the future; their distance will grow 
to about 8 x 10 6 km, approximately 20 lunar dis- 
tances, by the time of the first impact solutions, 
showing that its effect cannot be neglected in any 
future impact prediction. 
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5.3 Future Observations 

Our improved orbit for 2009 BD shows that it 
will pass at 0.9 lunar distances from the Earth on 
2011 June 2, becoming observable as it rapidly in- 
creases its solar elongation from 50° on May 28 to 
90° on June 1. Peak brightness occurs on 2011 June 
3 when the V magnitude reaches 17. 

Fig. 2 shows the ephemeris difference be- 
tween the gravity-only and gravity-plus-radiation- 
pressure orbit solutions during this observational 
opportunity. The peak difference of 9" coincides 
with the time of close approach, as the slowly 
increasing linear difference is modulated by the 
rapidly changing geocentric distance, creating a 
peak angular difference around perigee. The right 
ascension and declination components show differ- 
ent behavior with changing sign that should be easy 
to detect with careful astrometry. It should there- 
fore be possible to drastically improve the area to 
mass ratio determination. 

Other physical observations, including radar, can 
potentially provide constraints on geometric albedo 
and size, leading to a good density determination. 



Right Ascension 




A Peirce's Criterion (and its 
adaptation to a Rayleigh- 
distributed quantity) 

The following summary of Peirce's criterion 
is mostly based on the notation introduced in 
[9]; however, we ignore most of the discussion 
about the numerical iteration methods, that are 
not relevant nowadays in light of the current 
computational capabilities. The computation of 
the rejection coefficients is necessary in our case, 
since most precompiled numerical tables (from [9] 
or other sources) do not extend to the case of 7 
parameter fits, such as the one we need for our 
orbital computation (6 orbital elements plus the 
AMR). 

The basic idea of Peirce's criterion is similar 
to Chauvenet's: in a simplified way, we want to 
identify and reject any datapoint in our sample 
that lies more than a certain number of "error bars" 
(or sigmas) away from our fitted model, because 
we consider highly unlikely that such fluctuation 
may have happened because of random noise. 
The two methods differ only in how this "certain 
number", or rejection threshold, is computed* 15 ). 

In the basic version from Peirce, where the resid- 
uals are assumed normally distributed, the follow- 
ing relations allow us to compute the required rejec- 
tion threshold, via an iterative process: if we have 
N observations, of which n show doubtful residuals 
when fitted with an m parameter model, then the 
rejection threshold to reject all n observations is at 
x-sigma, where x is the positive solution of: 



Figure 2: Offset of the radiation pressure model 
from a purely gravitational solution, for 
the 2011 apparition. The plot starts on 
2011 May 27, at the time when the ob- 
ject will emerge from solar conjunction. 



x 2 = 1 + 



N — n — m 



Pi 
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R(x) n 
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' 15 'And for the already discussed fact that Chauvenet's is sta- 
tistically valid only for a single rejection, while Peirce's it- 
erative nature makes it usable even for multiple rejections. 
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where 



Pn,u = 



n 



N-n 



N -n 
N 



N 



and 



(p[x) 



R(x) = e 2 ip(x) 



e 2(f( = Erfc 
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As a side note, the rejection threshold for 
Chauvenet's criterion (let's call it x, to distinguish 
it from Peirce's x) could be obtained, in the 
above formalism, by simply solving the equation 

<p( X ) = Erfc (^) 



l 

2N 



With this formalism, a version of the criterion 
that is valid for our case of positive definite astro- 
metric residuals is easily obtained by replacing the 
Gaussian function in <p(x) with its Rayleigh coun- 
terpart: 



<p'{x) 



f 

J x 



2£e"« di = e 



It is interesting to note that, in this case, the 
function <p'(x) is analytically integrable, and it 
turns out to be a simple exponential function. 
This fact does not help the rest of the compu- 
tation of Peirce's threshold, that remains itera- 
tive and non-analytical, but it has interesting im- 
plications for Chauvenet's criterion: if we solve 
27v = f'ix) = /~ 2£e~ f2 (i£ = e~* 2 we obtain an 
analytical expression for the Chauvenet's threshold 
of the Rayleigh case, as 

B Debiasing and averaging 
method to extract the best 
value for H from sparse MPC 
photometry 

The following text explains in detail the pro- 
cedure used to extract a value for the absolute 
magnitude H of 2009 BD from the sparse and 
uncalibrated photometry made available together 
with the MPC astrometry. The proposed method 



(16) 



deals with debiasing, passband correction, rejec- 
tion and weighting, and it tries to make the best 
possible use of all the available data. 

The procedure starts with an average of all mag- 
nitude values of 2009 BD reported by each station 
in each night. This averaging is done first to reduce 
the noise in each datapoint, and second to properly 
take into account cases where a single station re- 
ported a large number of photometric points that 
might otherwise dominate the the final average and 
torque it toward a possibly biased value. This way, 
we obtain a total of 28 photometric datapoints, in- 
cluding 4 from our own observations from Mauna 
Kea. 

At the same time, the standard deviation of 
each set of points is also computed, and used 
as a first guess of the photometric error of the 
associated datapoint. This "statistical" error 
bar is then added in quadrature to a value of 
0.10 magnitudes; the reason for this correction is 
twofold: first, it takes into account a "quantization 
noise" introduced by the fact that each photo- 
metric point in the Minor Planet Center (MPC) 
format is given to a 0.1 magnitude precision, and 
second it places a lower limit on cases where the 
"statistical" error bar has a formal value of zero 
(because a single value was used in the average, or 
all values were equal within the reported precision) . 

Next we need to take into account the differ- 
ent passbands reported by each station, to con- 
vert them to a close approximation of the standard 
V magnitude, that is assumed by definition when 
computing the H value. At the same time, we need 
to compensate for known biases in the photometry 
associated with each station. Both these steps are 
accomplished during the same process, based on 
the historical performance of each site, as obtained 
from the AstDyS website^ 16 ). 

For each station that reported photometry of 
2009 BD, we downloaded from the above mentioned 
website the complete listing of all their observations 
of multi-opposition objects. The AstDyS summary 
file contains an estimate of the bias of each pho- 
tometric point, computed a-posteriori from the av- 
erage absolute magnitude of the object. These bi- 
ases take into account a conversion of each pass- 
band to V, obtained applying a predefined color 
thought to well represent the average color of an 

hamilton . dm . unipi . it/ astdys/ index . php?pc=2 . 
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asteroid in that passband. AstDyS uses a correc- 
tion V — R = +0.40 to correct for R- filtered magni- 
tudes, and the same correction is also used for the 
"placeholder" filter code C, that was allowed by the 
MPC at the time of the observations (and stands 
for "Clear", meaning no filter used)( 17 ). 

We restrict our analysis to those observations 
obtained in the years 2008 and 2009 (to account 
for possible changes in the instrumental setting of 
each observatory) , with the same passband and re- 
ported by the same program code that measured 
2009 BD (when more than one code is present at a 
station). We extracted the bias associated to each 
observations of these multi-opposition objects, av- 
eraged together those referred to the same object 
and the same night, and then averaged all the re- 
sulting values to get an estimate of the typical bias 
of the station. 

At the same time, we obtained an estimate of 
the noise on the bias, measured by the standard 
deviation of these single-object values. This 
scatter is mostly due to the rotational lightcurve 
of each asteroid included in the sample; if we 
estimate that there is an equal likelihood that each 
object was observed at any phase of its rotation, 
we may expect that these effects tend to cancel 
out in the average, decreasing their influence on 
the overall bias estimate as y/n^, where n s is the 
number of observations involved in the average 
of each station. Therefore, our best guess for the 
error bar associated to each station-specific bias is 
the standard deviation of its single biases, divided 
by y/n^; since the n s of observations involved in 
the average is quite high for most stations (more 
than 10 5 for the surveys, and around 10 2 for most 
amateur stations) we obtain error bars on the 
biases that are usually around 0.05 magnitudes or 
less. 

These values of the station-specific bias are 

'A correction of —0.80 is applied by AstDyS to the cases 
where a blank is reported in place of a filter code; this is 
based on the assumption (still in use at the MPC) that a 
blank filter code stands for a B passband (a standard in- 
herited from the photographic plate time). In most cases, 
observers reporting "blank" magnitudes nowadays mean that [ 
no filter was used, and therefore the assumption may be ex- 
tremely incorrect, introducing an artificial offset in the pho- 
tometry. However, the method used for our analysis will not 
be affected by this assumption, since this erroneous artificial 
offset will show up in the final average bias of the station, 
and will then be removed by the magnitudes we use in our 
average for H, that were also corrected under the same as- 
sumption. 



then removed from each photometric datapoint of 
2009 BD, and the two error bars (the one from the 
photometry, and the usually small one from the 
bias) are combined in quadrature, to get an esti- 
mate of the overall quality of the datapoint. These 
will be the starting values for our computation of 
H. 

The first pass of the H computation is a 
weighted fit of these corrected photometric data- 
points in the standard H-G formalism. Since the 
quality of the dataset in use is not very high, we 
are forced to assume a fixed value of G = 0.15, and 
restrict our analysis to a single free parameter H. 
This first fit, computed without rejections, gave 
a best-fitting value of H = 28.40. Residuals of 
each photometric point with respect to this value 
were then computed, and they formed the ba- 
sis for the rejection process, using Peirce's criterion. 

The application of the criterion to the magni- 
tude dataset is straightforward from a mathemat- 
ical point of view, since these residuals can be as- 
sumed to be normally distributed. However, a fur- 
ther physical complication is introduced by the fact 
that each residual is composed by two different 
parts, one given by the statistical noise, and one by 
the intrinsic lightcurve of the object, of unknown 
amplitude, period and phase. If we decide to reject 
a point based on the fact that its residual is a given 
number of sigmas away from zero, this may be due 
to the point being located at a maximum (or mini- 
mum) of the lightcurve, and not to its low quality. 
Because of the lightcurve, even a very high-quality 
photometric point may be rejected by a straight- 
forward application of the rejection criterion, if it 
lays at a lightcurve extreme. 

A first way to partially address this problem 
would be to artificially inflate the error bars of 
each data point, adding the lightcurve amplitude 
to each photometric error. However, this would re- 
quire an independent knowledge of the lightcurve 
amplitude, which is not available in this case^ 18 * 1 . 

For this reason, we have to use an alternative, 

Some uncalibrated photometry of 2009 BD, ob- 
tained from station 095, was made available online at 
www.astroalert.su/2009/01/18/neo2009bd; it shows that 
the object may have a lightcurve variation of approximately 
1 magnitude, and it seems to be variable on a timescale 
of 10 minutes or less. A clear evidence for lightcurve 
variation on the same timescale, without a period and an 
amplitude, can also be inferred from the data from station 
H10 mentioned in the text. 
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indirect method to extract the same information 
from the available data, in a statistical sense. If we 
assume that the residuals, after the first-pass fit of 
H, are due to the sum of the two effects (statistical 
error and lightcurve), we may expect that they 
are, on average, the sum (in quadrature) of these 
two components. We may then try to recover the 
unknown one of these components (the lightcurve 
amplitude, A;) from the other component (the 
statistical error, A s ) and the total error (measured 
by the residuals). If we do so on the average of all 
points, we see that the average statistical error is 
around A s ~ 0.16 magnitudes, while the residuals 
are scattered with RMS ~ 0.41 magnitudes. If we 
assume, m symbols, that A s + Aj 2 = RMS 2 , then 
we get Ai ~ 0.38, our estimate of the lightcurve 
amplitude^ 19 ). 

We can now add this estimate of the lightcurve 
amplitude to the error of each data point, and use 
this new error bar to estimate how many sigma 
each residual is away from the average. This is 
what we need to perform a rejection with Pierce's 
criterion, as discussed above. The test rejects 
n = 3 photometric datapoints, and the remaining 
25 values are then used to recompute a final value 
for H, and its error bar. 
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This lightcurve amplitude is estimated in an "RMS" sense, 
and does not correspond to the peak-to-peak amplitude usu- 
ally quoted for other objects. A conversion between the two 
may be obtained only under the assumption of a functional 
form for the lightcurve itself: if we consider a purely si- 
nusoidal case, then the peak-to-peak amplitude is equal to 
2\/2 ~ 2.83 times the RMS value, or about 1.07 magnitudes 
in our case, in good agreement with the estimate from sta- 
tion 095 mentioned above. 
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